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Abstract. The carrier dynamics in bulk Silicon, a paradigmatic indirect gap semiconductor, 
is studied by using the Baym-Kadanoff equations. Both the electron-electron (e-e) and 
electron-phonon (e-p) self-energies are calculated fully ab-initiohy using a semi-static GW 
approximation in the e-e case and a Fan self-energy in the e-p case. By using the generalized 
Baym-Kadanoff ansatz the two-time evolution is replaced by the only dynamics on the 
macroscopic time axis. The enormous numerical difficulties connected with a real-time 
simulation of realistic systems is overcomed by using a completed collision approximation that 
further simplifies the memory effects connected to the time evolution. The carrier dynamics 
is shown to reduce in such a way to have stringent connections to the well-known equilibrium 
electron-electron and electron-phonon self-energies. This link allows to use general arguments 
to motivate the relative balance between the e-e and e-p scattering channels on the basis of the 
carrier energies. 



1. Introduction 

Since the seminal works of two chemists, Norrish and Porter (recognized by the 1976 Nobel 
prize in Chemistry), laser technology has grown explosively until reaching in recent years laser 
durations of few attoseconds and power flux densities larger than 10^*^ W/cm^ [31j. This led to 
the opening of the new era of femtosecond (fs) nanoscale physics, where combined non-linear and 
non-adiabatic phenomena [46j occur and where it is possible to monitor real-time the dynamics 
of photo-induced electronic excitations with unprecedented precision [2T] . 

Ultrafast science at the nano-scale is clearing the path for nanostructure devices with efficient 
light emission spectra, high optical gain and controllable atomic deformations, thus paving the 
way to strategic applications in chemistry, biophysics and medicine. The design of nano-scale 
devices is, however, inevitably linked to a detailed knowledge of the structural and dynamical 
properties of these systems. Despite the massive number of available experimental results there 
are still scarce numerical and theoretical methods in use of the scientific community. 

Nanostructures and biological systems are, indeed, formed by hundreds/thousands of atoms 
and their peculiar properties are related to their reduced dimensionality and extended surface. 
Any reliable theory is inevitably linked to a detailed knowledge of their structural and dynamical 



properties. Due to the complexity of these systems, however, state-of-the art methods are 
confined to rather simple models with empirical parameters. In this way the theory is deprived 
of its predictive aspiration, and of the possibility of inspiring new experiments and practical 
applications. This situation is unavoidably creating a gap between theory and experiment. 

Density Functional Theory (DFT) [l5] represents the most up-to-date, systematic and 
predictive approach to study the electronic, structural and optical properties of an impressive 
wide range of materials ^30j, including solids and nanostructures. DFT is, however, a ground- 
state theory that, as it will be shortly discussed in this paper, can only provide a suitable, 
parameter-free (and thus predictive) basis where the actual real-time simulations are carried 
on. All the class of methods developed starting from the DFT are labeled as ab-initio methods. 

It is well known that extended systems have been mostly studied by using the 
many-body perturbation theory (MBPT) [35] and Time-Dependent Density-Functional- 
Theory (TDDFT) [39]. The different treatment of correlation and nonlinear effects marks 
the range of applicability of the two approaches. The real-time TDDFT makes possible to 
investigate nonlinear effects like second harmonic generation [52J or hyperpolarizabilities of 
molecular systems p2]. However the standard approaches used to approximate the exchange- 
correlation functional of TDDFT treat correlation effects only on a mean-field level. As a 
consequence, while finite systems — such as molecules — are well described, in the case of extended 
systems — such as periodic crystals and nano-structures — the real-time TDDFT does not capture 
the essential features of the optical absorption [l2] even qualitatively. On the contrary MBPT 
allows to include correlation effects using controllable and systematic approximations for the 
self-energy S, that is a one-particle operator non-local in space and time. 



2. The equilibrium and out of equilibrium dynamics in an ab initio framework 

In recent years, the MBPT approach has been merged with DFT by using the Kohn-Sham [39] 
Hamiltonian as zeroth-order term in the perturbative expansion of the interacting Green's 
functions. This approach is parameter free and completely ab-initio, |42] and in this paper 
will be addressed as ab-initio-MBPT (^i-MBPT) to mark the difference with the conventional 
MBPT. 

In practice the Kohn-Sham equation is first solved: 
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(t>nk (r) = Cnk'/'nk (r) , (1) 



with 



Vscf [{R}] (r) = V [{R}] (r) + J dr' v (r, r') n (r') + V^c ( [n] , r) . (2) 

Vs is, thus, written in terms of the bare Coulomb interaction u(ri,r2) = |ri— r2p"'^, the 

electronic density n(r), the exchange-correlation potential V^cli'^lji") and the total ionic 
potential V [{R}] (r): 

V [{R}] {r) = J2vir- R/) , (3) 

with R/ the position of the generic atom /. By plugging into Eql3]the precise atomic structure 
of a given material, the solution of Eq{T] provides a basis for the electronic levels suitable to 
rewrite the Many-Body problem. Indeed the lesser Green's functions can be easily rewritten as 
a matrix in the KS basis 

G< (1, 2;f3) = Y^ M iti,t2; /5)] Ck (ri) -/-n'k (ra) , (4) 

rm'k 



where 1 = (ri,ii) and 2 = (r2,t2) represent global position, time and (eventually) spin indexes. 

However the ^i-MBPT is a very cumbersome technique that, based on a perturbative concept, 
increases its level of complexity with the order of the expansion. As an example, this makes the 
extension of this approach beyond the linear response regime quite complex, though there have 
been recently some applications of the ^i-MBPT in nonlinear optics [13 \ \35 \ [36]. 

Another stringent restriction of the ^i-MBPT is that it cannot be applied when non- 
equilibrium phenomena take place: for example it cannot be applied to study the light emission 
after an ultrafast laser pulse excitation. A generalization of MBPT to non-equilibrium situations 
has been proposed by Kadanoff and Bavm [2^ Wl\ [5]. In their seminal works the authors 
derived a set of equations for the real-time Green's functions, the Kadanoff-Baym equations 
(KBE's), that provide the basic tools of the non-equilibrium Green's Function theory and allowed 
essential advances in non equilibrium statistical mechanics. 

In the equilibrium MBPT, due to the time translation invariance, the relevant variable used 
to calculate the Green's functions is the frequency u. Instead, out of equilibrium, the time 
variables acquire a special role and much more attention is devoted to the their propagation 
properties. The time propagation avoids the explosive dependence, beyond the linear response, 
of the MBPT on high order Green's functions. Moreover the KBE's are non-perturbative in the 
external field therefore weak and strong fields can be treated on the same footing. 

The key ingredient of the KBE's are the two times lesser /greater Green's functions [291 [25l 
MM G>„^k(*i'*2) and self-energies E>„^j^ (ti,t2)- If we consider the time evolution of 
only on the macroscopic axis (ti = t2 = T) we get the starting form of the KBE's I will further 
analyze in the later sections of this paper: 

i-^G< (T, T; /3) = [h^^ + V^^*™^ [n (T; - V^^*''- [no] + 

(T; 13) + Uk (T; /3) , G< (T, T; /?)] + S< (T; /3) . (5) 

In EqU] all quantities are matrices in the space of KS states (see EqiH), e.g., S^^,^^,^, /3 = 

(fc^Tg;)^^ where T^i is the phonon bath temperature, and I have already introduced the 
fundamental ingredients that will be discussed in the following of this paper: 

• The KS Hamiltonian, h^'^. 

• The Hartree potential V^"''*''^^ The ionic and kinetic part is already embodied in h^'^. 

• The static self-energy S^. This is responsible for introducing static e-e correlations that, 
if properly defined [Ij reduce the KBE, in the linear-response regime, to the well-known 
Bethe-Salpeter equationjM]. 

• The interaction with the external, arbitrary electric field. The only assumption here is the 
the field wavelength is long compared with the unit cell size (optical limit). 

• The interaction kernel S^. This is the most important part of the equation that embodies 
dynamical e-e and e-p scatterings. It will be defined and carefully analyzed in SecjH 

One of the first attempts to apply the KBE's for investigating optical properties of 
semiconductors was presented in the seminal paper of Schmitt-Rink and co-workers. [15] Later 
the KBE's were applied to study quantum wells, [43j laser excited semiconductors, [2B] and 
luminescence [22]. However, only recently it was possible to simulate the Kadanoff-Baym 
dynamics in real-time. [Ml SH (Uj 

Regaring the out-of-equilibrium carrier dynamics there is a vast litterature [42| where the 
KBE's have been solved using several kind of approximations. I do not want to give here an 
exhaustive list of all the works in this field, but a key aspect of most of them is that they rely on 
very strong approximations for what regards the underlying electronic structure and microscopic 



dielectric properties of the materials used. For example a large number of works use few (two 
or four) bands models like in Refs.|17l 1511 \53\ \28[ [6l [15]. In all these works the electronic levels 
are approximated as well as the screening of the electron-electron interaction. 

In this work I want to go well beyond these approximations by solving the KBE's in a fully 
ab-initio iiamework. I will refer to this novel methodology as a6-miiio Non-Equilibrium Green's 
function (^i-NEGF) theory. 

3. The equilibrium self energies 

In order to discuss the physics involved in the dynamics of carriers taken out-of-equilibrium it 
is useful to introduce the properties of the two self-energies that will be used to describe the 
e-e and e-p coupling. In the next two subsections, indeed, I will introduce the approximations 
used and the key definitions of the different ingredients needed to build the self-energies. I will 
later rewrite these equilibrium self-energies in the framework of the ^i-NEGF theory. 

The e-e interaction will be treated in the so called GW approximation [5J where the self-energy 
is expanded at the first order in the dynamically screened interaction. This approximation 
has provided a priceless tool to evaluate accurately the electronic levels of a wide class of 
materials |42] . 

Since its first application to semiconductors [50] the GW self-energy has been shown to 
correctly reproduce quasi-particle energies and lifetimes for a wide range of materials [5j. 
Furthermore, by using the static limit of the GW self-energy as scattering potential of the Bethe- 
Salpeter equation (BSE) [39], it is possible to calculate response functions including electron-hole 
interaction effects. 

In contrast to the e-e case the e-p coupling has received only minor attention in the ab- 
initio community. And this is somehow astonishing if we consider that the e-p coupling is well 
known to play a key role in several physical phenomena. For example, beside inducing the 
well-known superconducting phase of several materials [H], it affects the renormalization of the 
electronic bands [IJ [IH [Ml [El [TO] 

Despite the development of more powerful and efficient computational resources the actual 
calculation of the effects induced by the e-p coupling in realistic materials remains a challenging 
task. In addition to the numerical difficulties, it has been assumed, for a long time, that 
this interaction can yield only minor corrections (of the order of meV) to the electronic 
levels. As a consequence the majority of the a6-imiio simulations of the electronic and optical 
properties of a wide class of materials are generally performed by keeping the atoms frozen in 
their crystallographic positions. More importantly it is well-known that phonons are atomic 
vibrations and, as a such, can be easily populated by increasing the temperature. This naive 
observation is de-facto used to associate the effect of the EP coupling to a temperature effect 
that vanishes as the temperature goes to zero. However this is not correct as the atoms posses 
an intrinsic spatial indetermination due to their quantum nature, that is independent on the 
temperature. These quantistic oscillations are taken in account by the e-p coupling when T — )• 
in the form of a zero-point-motion effect. 

3.1. The GW approximation for the electronic self-energy 

The Feynman diagram corresponding to the GW approximation (GWA) is showed in the panel 
(6) of Fig[Tl It corresponds to the lowest order contribution to the self-energy expressed as a 
Taylor expansion in the screened interaction W (1, 2): 



In the present case we will use a spinless formulation where spin indexes are summed out. The 
extension of the full theory (at the equilibrium and out-of-the-equilibrium) is straightforward. 
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Figure 1. The self-energy diagrams used to account for e-e and e-p coupling. For the e-p case I consider the 
first order expansion in the atomic displacements corresponding to the Fan self-energy. In the e-e case, instead, 
I first consider an equilibrium GW self-energy (diagram (b)) and the corresponding out-of-equilibrium version 
(diagram (c)). In this last diagram the independent particle polarization corresponds to the dynamical part of 
the diagram. Indeed the two wiggled lines correspond to a statically screened microscopical e-e interaction. 
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Figure 2. The equation for the screened interaction W within the RPA. 



In the GWA the screened interaction is written as: 



W (1,2) = V (1,2) + I d34t;(l,3)x(3,4)VF(4,2) 



(7) 



with f (1,2) = V {ri,r2) 6 (ti — is the bare Coulomb interaction and x is the irreducible 
electronic response function [l9]. An approximation for W follows by choosing x to be calculated 
in the independent-particle approximation which corresponds, for W, to the random-phase- 
approximation (RPA) , where 



X (1, 2) « xo (1, 2) = -iGo (1, 2) Go (2, 1) . 



(8) 



The actual evaluation of the space dependence appearing in EqiH] is treated in reciprocal space, 
by summing on the lattice vecors G, G' and integrating over the Brillouin Zone (BZ): 



X(ri,r2;t) 



BZ (27r)^ 



(9) 



It follows that, within the Random-Phase approximation (RPA) : 



xg g' (q, ^) = xg g' (q. ^) + ^xg g" (q. 
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EqJlOlis used together with the definition of the single-particle Green's function 
Go (ri,r2;a;) = 2^ ^ </>„k (ri) Ck (^2) G„k (w) , 



(10) 



(11) 



with 
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(12) 



to write the final expression of the GW self-energy . In the present work we are mainly interested 
in the carrier dynamics and, in order to keep the notation as compact as possible, we restrict our 
attention on the imaginary part only of Sg^^g. This, indeed, defines the quasi-particle linewidth 

as 
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where we have used the simplest on-the-mass-shell approximation [5]. By doing simple 
algebraical operations it is possible to rewrite as: 



I'nk"'''^ = ^ X] X] X] ^' [^"™ ('^' ^')] * [^G,G' (q, enk - Cmk-q)] 



m q G,G' 



(2 /»nk— q) ^ (^nk ^mk- q) /mk— q^ (^mk- q ^nk) 



electron scattering 



hole scattering 



■ (14) 



The term (electron scattering) represents an electron with energy e„k that decays into an 
empty state with lower energy e^k-q dissipating the corresponding energy difference in positive 
energy e-h pairs (represented by e~^). The {hole scattering) channel is equivalent to the 
{electron scattering) with the only difference that the hole jumps into a full state with higher 
energy. 

In EqUHwe have defined 



(k,q,G) = (nk|e*('i+*^) ''i|mk-q) = J drul^ (r) 



Umk-a (r) e 



i{q+G)-r 



(15) 



where n„k (r) is the periodic part of the Bloch function. We have also introduced the dynamically 
screened interaction VFg,g' (q^^^) defined as 



w^G,G' (q,w) 



(47r) ^ (47r) , , (47r 
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(16) 



3.2. The electron-phonon Fan self-energy 

The phonon contribution to the equihbrium self-energy comes from the oscillations of the atoms 
around their equilibrium positions. By starting from Eq|T] we divide the total Hamiltonian of 
the system in electronic {Hei), atomic {Hat) and electron-atom part (Hei-at)'- 

H = Hel + Hat + Hel-at, (17) 

with 

Hei-at= [ drp{r)Vscf[m]{r). (18) 

J crystal 

The vibrational states of the Hamiltonian H are described, fully ab-initio , by using the well- 
known extension of DFT, Density Functional Perturbation Theory (DFPT) [71 [20] where the 
change in the ionic potential due to a generic atomic displacement is taken in account, self- 
consistently, with the solution of the KS equation. Indeed the total electronic potential depends 
on the atomic positions R/. 

If we now consider a configuration of lattice displacements u/^, H can be expressed as a 
Taylor expansion 



H-H-l^ MTa 

la 



where a is the Cartesian coordinate. Any overlined operator O is evaluated with the atoms 
frozen in their crystallographic equilibrium positions. The link with the perturbative expansion 
is readily done by transforming Eq. ()19p from the space of the lattice displacements to the space 
of the canonical lattice vibrations by means of the identity [10]: 

iiia = i2N,Mju;^xr'^' (qA|/) e^^^ ^^^) (bl^^ + b^x) , (20) 

qA 

where is the number of cells (or, equivalently the number of q-points) used in the simulation 
and Mj is the atomic mass of the I-th atom in the unit cell. i^M^) is the phonon polarization 
vector and Sl^^^ and 6qA are the bosonic creation and annihilation operators. 



By inserting Eq. (j20l) into Eq. ([T9|l we get 

H-6=^ ^^'-Ik^Wnk p.q, + Q (21) 



/Nq 

krm qA 



In Eq. ()2ip we have introduced the (s^/^j^) electron-phonon matrix elements. This is obtained 
by rewriting Vscf and by making explicit its dependence on the atomic positions: 



t>e/[{R}](r) = j;v;g(r-R,,). (22) 



Is 

By using Eql22II get an explicit form for the electron-phonon matrix elements 



QV^''^ (r) 

= E m.u^xr'^^ («'k + q| \nk)o^a (qA|s) , (23) 



with the r integral appearing in (n'k + q| — |?^k)o performed in the unit-cell instead that 
in the whole crystal. We have used the short form Rsa = Risa\i=Q- 

The actual calculation of the Fan diagram is straightforward. The Fan diagram is similar to 
the GW one where the screened electronic interaction W is replaced by a phonon propagator 
of wavevector q and branch A|37j. By applying the finite temperature diagrammatic rules it 
is possible to define the Fan self-energy operator Yf^{uj)^ recovering the expression originally 
evaluated by Fan [16] for the phonon mediated quasi-particle linewidth 
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with Nf^x (T) is the Bose function distribution of the phonon mode (q, A) at temperature T. 
The e-p lifetime is splitted in Eg 125) in three different contributions corresponding to the 



scattering of an electron/hole F^j^ 



e-p,eq 



e/h 



with the phonon bath and an explictly temperature 



term, F^j^^'^^ (/3) . The Fan self-energy can be also derived by using the framework of the 

ph 

Heine- Allen-Cardona (HAC) theory [21 [U |Tl]. The HAC approach is based on the static 
Rayleigh-Schrodinger perturbation theory done using the ujsa displacement operators as scalar 
variables. 



3.3. Electron-phonon versus electron-electron contribution to the equilibrium electronic 
lifetimes 

From Eq(T3]and Eql25]we can deduce a simple physical property of the equilibrium self-energies. 
By simple energetic conditions it is clear that in a system with a direct or indirect gap Eg the 
e-e F^^^'*^'' are exactly zero whenever Cnk < ^cbm + 2£'g or e„k > ^vbm — 2-E'g. These two 
energy conditions define the energy ranges represented by gray areas in FigiSj 

On the contrary the e-p linewidths, thanks to the continuum of phonon frequencies, have no 
zero energy regions. It follows that in the case of low energy electrons or hole the e-p linewidths 
largely dominate on the e-e part. Only above (or below) the energy forbidden regions the e-e 
linewidths start growing quadratically, following the Fermi liquid expected behavior, and cross 
the e-p part at about — 4eV and 6eV. 



4. An out of equilibrium formulation of the carrier dynamics including e e and 
e p scatterings 

Now that the equilibrium self-energies have been introduced and defined we can move now in 
the out-of-equilibrium regime by giving a precise form to the interaction kernel Sk {T) . First of 
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Figure 3. Quasiparticle linewidths of bulk Silicon calculated by using the GW approximation for the e-e 
scattering (green boxes) and the Fan approximation for the e-p scattering (red circles). The two gray areas denote 
the energy regions (as large as twice the electronic indirect gap 2Eg «1.2 eV) where, by simple energy conservation 
arguments, the electronic linewidth are zero by definition. In these energy conditions the e-p contribution is 
stronger and the corresponding linewdiths are larger then the e-e ones. The quadratic energy dependence of the 
e-e linewidths inverts this trend at about five times the indirect gap. 

all this is formally defined [291 [251 iZl [H] as 

Sk (T; /3) = r dr (T, r; /3) G< (r, T; /3) + G< (T, r; f3) S> (r, T; f3) - 

J —OD 

S< (T, r; /3) G> (r, T; /3) - G> (T, r; /3) I]< (r, T; /3)] . (26) 

The key point in Eq[26|is that the definition of S is not closed as it involves two-times function 
and both the greater and lesser Green's functions. This means that one should consider the two- 
time evolution of the KBE together with the equation of motion for the G^ Green's function 
and, as we will shortly see, G^ the corresponding retarded quantities. 

As a two-time evolution of the KBE is, at the present stage, computationally not feseable 
for the real-time dynamics of complex solids and nanostructures we want to explore a different 
path based on several key ansatz and approximations. 

For the e-e case the extension of the GWA, EqjTl to the out-of-equilibrium regime is well- 
known |25j. The diagrammatic form is represented in the frame (b) of Fig[Tl What I want to do 
here is the rewrite this extension in a a6~inzizo framework introducing some key approximations. 
First of all I split S and S in the e-e and e-p parts: 

ih,t2;P) = ^t"^^ ih,t2; /3) + (ti, t2; /3) , (27) 



(T; /3) = Sj,'-'^^ (T; /?) + S|,'-^^^ (T; (3) . (28) 



The GWA for out-of-equilibrium reads: 



nil ,7712 Gi,G2 q 

Wla, (q, ti^t2;f3) *2; /3)}, (29) 

with the screened interaction written as 



g;3,G4 



^G3,G4 ^1' /3) T^g,G2 -^2, t2; /?) . (30) 

In EqlSO] Vl^(''"/") are the retarded/advanced screened interactions and is the independent 
particle polarization diagram represented by the closed loop in the diagram (b) of FigH} 

^Gi,G2 (q'*l'*2;/3) = Y Y [Pmam'a (P,q,Gi) Pm>2(P'q'G2) 

, {ti,t2\ 0) , (t2,ti;/3). (31) 

Before analyzing in detail the temporal dependence of it is instructive to introduce the 

out~of-equilibrium e-p self-energy. The extension out-of-equilibrium of the FAN self-energy is 
straightforward [25] : 

^Tn'f (ii>i2;/3) = 9 Y E(<'n.ik)%™k^!(q'*i'*2;/3)G^,^,k-q(ii>i2;/3)}, (32) 

mi,rra2 q,A 

with Dx the phonon propagator that I will shortly define. 

Now I introduce the complete set of approximations I will use to transform Eql2S]in a simpler 
set of equations of motion to describe, ab-initio , the carrier dynamics in realistic materials. 
These approximation are guided by simple physical arguments introduced with the aim of 
limiting the computation effort. This is a crucial requirement to make the simulations feasible. 
At the same time those approximations represent a first step in the device of a sound and 
complete framework for the ^i-NEGF approach and their validity will be carefully analyzed in 
the further development of this field. 

The approximations are: 

(i) The screened interaction entering Eql30] are assumed to be static. I introduce this semi- 
static approximation (SSA) of the out-of-equilibrium version of the GWA being inspired 
by the well-known Coulomb-Hole plus Screened- Exchange approximation (COHSEX) [52]. 
The SSA reads: 

^g/,G2 (q' ^1 ' /?) « ^G,,G2 (q) • (33) 

This approximation is well motivated in semiconductors and insulators as an electronic gap 
of 1 eV corresponds to a time scale of 0.66 fs. This a very short time compared to the time 
scale of the simulations presented in the present work. 



(ii) The simulation is done at a fixed temperature (3 that is assumed to not change during 
the simulation. This approximation corresponds to take the phonon population at the 
equilibrium and assuming that their energy distribution does never deviate from a Bose 
distribution. In practice this means that: 

D< (q,ti,t2;/3) « (-^)E<A (/3)e±-<^^(*^-*^). (34) 

± 

Note that D> (q, 1 1 , ^2 ; /3) = - [Z?^ (q, ti , t2 ; ^S)] * . 

(iii) I use the Generalized Baym-Kadanoff ansatz (GBKA) |25 | H7 1 [H]. 

(iv) I neglect the off-diagonal matrix elements of the and 5]^ matrices 

{T,T-P)^ ibnn' /nk (T; /3) , (35) 
Sfn'k /5) ~ -^nn'Snk (T; /3) . ([33) 

As it will be showed in the following these approximations will make possible to perform the 
actual simulations of complex materials at a reasonable computational cost. 



^.i. The collision operator in the Generalized Baym-Kadanoff ansatz 

The GBKA provides the most important tool to disentangle the complex two-times evolution 
by rewriting: 



G^k *2; /3) « (^) (ti - /3) (*2, t2;f3)- (^i^ /3) G'2 (ti - ts; /3) 



(a) 



(36) 



that, in the lesser case, acquire a simple form written in terms of the time-dependent electronic 
occupations: 



G'<k iti,t2;P) « - g2 (ti - 12; /3) /nk (t2; /3) - /„k (ti; /3) G^^i^ {h - ts; /?) 



(a) 



(37) 



A similar identity holds for the greater Green's function. By using EqJ36] the KBE's for 
the electronic occupations /„k {T, (3) acquire a closed form if a given analytic form for the 
advanced/retarded Green's functions is given. This is the present case. But before discussing 
the specific form of the (r/a) Green's functions I use Eqi36]to rewrite the two contributions to 
S. 

The algebraical passages to rewrite S^^^'^ by using the GBKA can be found, for example, in 
Ref . |23[ [M] in the simpler case of a Jellium model with e-p interaction treated using a Frohlich 

interaction ff^^^ ~ fl'lql ^^"^ ^ single optical phonon mode. 

After some involved but conceptually simple algebra it can be shown that 
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As it will be clear in the next section the sum of S^^'^ and of its adjoin in Eg 1551 is crucial to 
keep the time-evolution unitary in such a way that the number of electrons is conserved. The 
actual time decay of the photo-excited carriers will, then, be connected to the actual memory of 
the system linked with the r integral appearing in Eql38l rather that with the imaginary part of 
the self-energy, as commonly done in equilibrium theories. This detail marks the basic physical 
difference between the equilibrium QP linewidths and the out-of-equilibrium QP decay time, as 
it will be clear from the discussion in SeciSl 

In the case of the e-e scattering the application of the GBKA is more involved as it must 
be applied to three Green's functions instead of one, like in the e-p case. Nevertheless, once a 
simple form for G^''/") will be given the physical interpretation of the collision operators will be 
equally simple. By following a procedure similar to the e-p case we find that 
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Eq J38fH0] represent the expressione for the time-dependent kernel S that I will further simplify 
in the next section. By using this expression for S the final equation of motion for the electronic 
occupation reads: 
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The Completed Collision Approximation 
Thanks to the GKBA the BKE's can be rewritten as a set of non-linear differential equations 
for the occupation functions. Nevertheless the whole history of these functions should be 
integrated numerically in Eqs l39lH0l In order to further simplify the theory I introduce two 
more approximations to the list presented in Sec 121 

(v) The functions are assumed to have a QP form (exponential): 
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With Gi^(r;/3) 



g2(-T;/3) 



This means that we neglect any deviation from the 

pure exponential behavior induced by the short-time dynamics [Q] and the e-e and e-p 
correlations [,24J. In eq|l2]the QP energies Enk and widths r„k are calculated ab-initiohy 
means of the equilibrium GWA and Fan self-energies. The explicit expressions for P^k are 
given by EqUH and Eql25l 



(vi) I adopt the so called Completed Collision Approximation (CCA) [23j that makes possible 
to solve the time-integration appearing in Eqs l39lH0] analytically by assuming that the 
occupation functions vary only slowly on the time scale defined by G^^^ (T — r; /3). This 

time scale is given by r„k (/3)~^ and, as shown in FigjSl it is of the order of 10 20 fs. Under 

this assumption the / functions can be taken out of the integrals appearing in Eqs l39lH0l 
and the rest of the integral can be done analytically by using Eg 1421 

Thanks to the CCA Eql39]can be rewritten as: 
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By using EqJ42lthe r integral reads: 
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Thanks to EqiMl the S^^' (T; /3) functions can be rewritten in terms of three simple lifetimes: 
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In order to define and interpret the physical properties of these lifetimes I introduce two 
functions: P"^* and 
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These two functions represent the energy constrain in the emission or absorption of a phonon 
smeared out due to the QP approximation for the retarded Green's function that. This, indeed, 
yields a Lorentzian indetermination in the energy conservation. By using these two simple 
functions the three e-p induced lifetimes can be easily rewritten as 
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The dynamics induced by the e~p kernel on the electronic occupations can be easily visualized 
using Eg J45I plugged inside Eql5] after taken the diagonal of . If I do not consider the coherent 
part of the KBE equation I get 
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EqH9] makes the interpretation of the carrier dynamics induced by the e-p coupling 
straightforward. We start indeed noticing that, simply because of their definition, 7^^^ ^'^^ and 

7nk ^'^^ always positive. In contrast 7^^ ^'^'^^ can be positive or negative. 

To see why this property is important let's consider two limiting cases when f^k (/3; T = 0) = 
or fnk {f3;T = 0) = 1. In the first case, from EqUO] it follows that, if we consider the zero 
temperature case. 
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(50) 

m) 



This simply means that an initially completely filled level will be emptied by the e~p interaction 
while an initially empty level will be filled. The 7^^ ^'^^^ is non zero only at finite temperature 
and works to build up a finite temperature quasi-equilibrium among the carrier occupations. 
The 7^^ P'^/'') lifetimes represent a single scattering of the state |nk) to the state |nk — q) with 
emission/absorption of a phonon mode |qA). This process is similar to the one studied in the 
equilibrium case and, indeed, if we consider the Lorentzian function appearing in P"^^ and p^"*** 
as approximated delta functions we see that the analytic form of 7^^ P''^/'*) g^j-g proportional to 
two terms appearing in Eql25[ Nevertheless, as we will clearly see in the case of bulk silicon, the 
dynamics induced by the KBE equation is a cascade process where many phonons are emitted 
or adsorbed and the final decay time of the electronic levels will be very different from the 
equilibrium case. 

By using approximations (vi) and (vii) the e-e kernel, EqUU] can be worked out in a way 
similar to the e-p case. Although the math is more involved the final outcome is quite similar 
to the e-p case. Indeed it turns out that 
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The two additional lifetimes 
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We have introduce the function R defined as: 
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The physical interpretation of the e-e assisted carrier dynamics is very similar to the e-p case. 
The only difference is that, instead of phonons, the carrier decay is due to the creation of e~h 
pairs. The basic transition that contributes to the e-e scattering is a four bodies process: the 
initial electron (or hole) scatters in an e-h pair plus another electron (or hole). If we distinguish 
between the bath electrons and the optically pumped ones (the carriers) we see that these four 
bodies can, quite generally, be distributed in both group of states. 

Nevertheless in the case of semiconductors and insulators it is meaningful to distinguish 
between intra-carrier and extra-carrier channels. In the first case the e-h pair is created among 
the optically pumped carriers while in the second case the pair is created across the gap involving 
the bath electrons. Of course these two channels, in general, contribute at the same time to 
the overall dynamics. But, as we will see in SeciS] their relative strength depends on the carrier 
energy. 

In conclusion the total rate equation governing the carrier dynamics that will be used in the 
Seci5]is given by EqHUwith defined by EqsHSlandlSTl I coded this equation in the Yambo 
code [3], in the same numerical framework where the equilibrium GW and Fan self-energies 
are coded. The existing structure of Yambo has already produced several publications and it 
provides, therefore, the most appropriate environment to perform out-of-equilibrium simulations 
in a controllable manner. The e-p matrix elements and phonon modes are calculated, instead, 
with the Quantum Espresso package [T8]. 

4-3. A detailed balance condition to ensure that the number of electrons is conserved 
In the case of e-e scattering the actual possibility of creating high energy electron-hole (e-h) 
pairs can modify, during the real-time dynamics, the number of carriers that fill the conduction 
bands with respect to the pairs initially photo-excited. This process is particularly evident for 
carriers created with energy larger than 2Eg below the VBM, or 2Eg above the CBM. 

In the e-p case, instead, the maximum available energy that any carrier can dissipate during 
the relaxation is the Debye energy that, in the systems I consider in this work, is much smaller 
than the gap. In this case, indeed, it is important that the theory preserves the correct number 
of photo-excited electrons during the whole simulation. 

I use this example to discuss a methodological aspect embodied in the use of a fully ab- 
initio iiamewoTk. Indeed the e-p matrix elements defined in Eq j23t can be written as 

9^L ^ 9n'nt = E (2M,a;k-pA)-^/' e^^'^-^y^' {n'p\^^\nkU^ (k - p, A|.) . (54) 

We can, now, consider the KBE for the two isolate states |nk) and |mp) that are scattered 
one into the other by means of the •S'^j^^''" For simplicity I only consider the dynamics 

induced by the e-p kernel at zero temperature and for a single phonon band ^. This dynamics 
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Now we notice that, thanks to simple symmetry arguments it can be shown that 
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Moreover from Eqs H6H171 it follows that -P^^**k_p)^ — ^mnp{p-k)i^i ^-nd, consequently, from 
Eqsl55l we see that: 
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By summing Eq 1571 and Eqi58]one gets 
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But from the boundary condition (2 — fmp {T = 0) — fnk (T = 0)) = it follows 
^ (2 — /mp (T) — /„k (r)) = 0. The total number of electrons in the two levels made scat- 
tering by the e-p interaction is conserved, as it should be. 

Eql56] thus represents a sort of detailed balance condition that ensures that the fraction of 
charge lost by a given state |nk) is exactly compensated by the charge added to the |mp) 
state. In total there is no unphysical charge lost or increase. Of course the level of accuracy 
of this detailed balance condition is ensured by the numerical precision under which Eqi56] is 
satisfied. This precision must be carefully checked and, eventually imposed, as on long-time 
simulations (of the order of picoseconds) tiny deviations can induce important breakdown of the 
total number of electrons conservation. 



5. An application in a paradigmatic solid: bulk silicon 

After all the math and approximations introduced in the previous sections it is instructive to 
apply all the machinery composed by the KBE's written in the ab-initio iramewoik to perform 
realistic simulations. The paradigmatic case I consider is bulk silicon. This a quite simple 
material that, neverthele^^ Inns an intprpstincr r»rnr»prt'\7"' tVip miTiimnl cran ig indirect. 




Figure 4. Time dependent occupations (lower frames) /„k (T; /3) as a function of the T for two external 
fixed temperature (0 K left frame and 300 K right frame) for three selected states in the band-structure of Silicon 
(upper frame). At T = only the Fis is artificially populated with 4 electrons needed to fill 2 of the 3 energy 
degenerated states. This explains why the y-axis goes from to 4. In the lower frames the dynamics is done 
including only e-p scattering (continuous lines) and e-e plus e-p scattering (dashed lines). In this case, as shown 
in Fig[31 the energy of the IFis) state is below the 2Eg threshold. This means, indeed, that at T = the e-p 
scattering is largely dominant over the e-p one. But, as the time increases and the initial charge is fragmented in 
lower energy states the intra-carriers e-e scattering quickly increases leading to a much faster population of the 
Xi state. 



From the point of view of the carrier dynamics this is an important aspect. Indeed the optical 
gap is located at the F point and it is larger (~ 3.2 eV) of the indirect gap (~ 1.2 eV). Thus 
if an electron is optically pumped in the conduction states this will populate the energy region 
surrounding the CBM, that is the F15 state, and the relaxation is expected to occur along the 
FX line (see the upper frame of FigS]) towards the Xi point. 

This has been, indeed, experimentally measured [27] by using Time-Resolved Two-Photon 
Photoemission (2PPE) and the time-decay of a level measured along the FX has been found to 
be around 40 fs. 

I do not want, however, to interpret and reproduce the 2PPE experiment but to use bulk 
silicon to illustrate the relative balance of the e-e and e-p channels in the carrier dynamics on 
the basis of the initial carrier energy and density. I will also briefly discuss the effect of the 



scattering with acoustic phonons keeping in mind the large number of simulations done using 
simple optical-only phonon models. 

The simulation is performed by defining an ad-hoc initial condition for the occupations and 
let the system evolve by means of Eql5j In the first series of calculations I take 

/nk (t = 0; /3) = 2 when |nk) = [Fis). (60) 

Two electrons (to account for the two spin channels) are moved by hand in two out of the 
three states that are energetically degenerate with the Fis stat. This corresponds to an initial 
carried density of 3.2 x 10^^ cm~^ . As the Fis state is at the border of the 2Eg threshold its 
extra-carriers contribution to the e~e lifetime is very small. This is because the energy of the 
state is not enough to create e-h pairs across the minimum gap. The real-time occupations of 
three paradigmatic states: jTis), \Li) and \Xi) are showed in FigUl In the k-grid I have used 
the Xi state is the nearest to the CBM. 

In FigjUthe functions /r^s {T;/3) (red curves), J'l^ {T;I3) (blue curves) and fxi {T;P) (green 
curves) are showed comparing the dynamics with only the e-p scattering (continuous lines) and 
with the e-p and the e-e scatterings (dashed lines) included. We immediately see some clear 
trends. The decay if the Fisis very fast, with an exponential time-dependence with a lifetime 
of the order of 12 fs. This value is very similar to the value of 16 eV obtained by calculating 
Fp~g '^'^ (/3 — )■ oo) showing that the e-e contribution is practically zero. As the time increases 
the initial charge is fragmented and lower state are populated. At this point the intra-carriers 
contribution to the e-e is rapidly increased and the contribution to the filling of the Xiand 11 
states cannot be neglected. Indeed by considering only the e-p channel the rapid increase in the 
occupations of the two states is largely underestimated. 

The message is thus clear. In the very short time the dynamics of carriers with energy below 
the 2Eg threshold is e-p dominated. Nevertheless in this regime the use of the CCA is highly 
questionable and a more careful full-memory treatment is required. As time increases also in 
this energy range the intra-carriers contribution to the carrier dynamics makes the e-e scattering 




Figure 5. Time dependent occupations (right frame) /„k (T; /3) as a function of the T for an external 
vanishing temperature for three selected states in the band-structure of Silicon (left frame). At T = only the 
Fis is artificially populated with 4 electrons needed to fill 2 of the 3 energy degenerated states. In the right frame 
the dynamics is done only including the e-p scattering but including only optical phonon modes (continuous lines) 
or all modes (dashed lines). We see that the effect of the acoustic modes is not at all negligible, especially for the 
Xi state filling. 

An important aspect that must be taken in account in the analysis of the time-dependencies 
showed in FigS] is that the overall filling of the Xiand 11 states is enhanced when the e-e 



scattering is included also because of the promotion of additional electrons from the valence to 
the conduction by the creation of e-p pairs across the gap. This process, although energetically 
impossible for the Fisstate, is made possible by the energy indetermination caused by the use 
of the QP approximation for the retarded Green's function. Indeed thanks to the Lorentzian 
line-shape, the R functions (see Eq l53p can be non-zero for arbitrary values of the energy levels 
involved in the e-e scattering. This is, in my opinion, a stringent limitation imposed by the QP 
approximation that can be easily avoided by using more refined approximations or by calculating 
the equilibrium spectral functions numerically in the GW and Fan approximations. 

By looking at the two lower frames of FigH] we see that the increase of the temperature 
from to 300 K does not appreciably change the occupation dynamics in the first 500 fs. This 
is reasonable as the e-e lifetimes are not explicitly temperature sensitive and, as far as the 
e-p channel is considered, the Silicon Debye temperature is as large as 645 K. Thus at room 
temperature only the low-energy acoustic modes are populated. Although those contribute for 
half of the e-p lifetime the overall effect on the dynamics is marginal. 
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Figure 6. Time dependent occupation as a function of the T for the set of electronic states with energy 
5.76 eV, well above the 2Eg threshold energy. The simulation is performed with the e-p scattering and including 
(dashed line) and not including (continuous line) the e-e contribution. For this high energy level the equilibrium 
e-e and e-p lifetimes are similar and both contribute to the decay of the state. This contrasts the Fiscase where 
the e-p scattering dominates in the short-time dynamics. 

An important issue that can be fully exploited by the present calculations is, indeed, the role 
played by the acoustic phonon modes that, in model calculations, are often neglected. As shown 
in Figj5] the acoustic modes, even at zero temperature, contributes for more than half of the 
filling of the Xi state. Their contribution to the Fisis, instead, minor and it slowly increase as 
the energy of the state approaches the CBM. This phenomena is connected with the physics of 
both the e-e and e-p induced relaxation of the carriers. This is indeed a cascade process where 
the initial Fiselectrons jump towards other states by the emission of, initially, a single phonon 
and a single e-h pair. This elementary process is very fast, of the order of even faster of 10 fs. 
But in order to reach the Xi state the electron undergoes several of these elementary processes 
and the overall relaxation is a much slower process. 

The last aspect I want to analyze by using bulk Silicon as paradigmatic system is the relative 
influence of the e-e and e-p scattering channels in the short-time regime for carriers pumped 
well above the 2Eg threshold. To this end I change the boundary condition give by Eql60]to be 



/nk (t = 0; /3) = 2 when E„k = 5.76 eV 



(61) 



This corresponds to an initial carried density of 1.9 x 10^^ cm~^, an order of magnitude larger 
then in the Fiscase. Prom Figl3]we see that the equilibrium QP e-e and e-p lifetimes are very 
similar. And, indeed, the decay of the 5.76 eV state is very much affected by the e-e scattering 
that, in this case, at short times is dominated by the extra-carriers channel. 

6. Conclusions 

In conclusion I have presented an ab-initio approach to the solution of the Baym-Kadanoff 
equations in realistic materials. The present approach is based on two major approximations: 
the generalized Baym-Kadanoff ansatz and the completed collision approximation. These makes 
the calculation feasible, as shown in the case of bulk silicon and provide a consistent and solid 
theoretical and numerical framework that can be applied to any kind of material that can be 
studied using the standard a6-inziio Many-Body Perturbation Theory techniques. Indeed, a 
major point of the present approach is that all the presented theoretical and numerical tools has 
been coded in an existing long-term, stable and highly productive code: the Yambo code[3|. By 
using bulk silicon as a paradigmatic material I show that, in semiconductors and insulators, the 
presence of an electronic gap makes the e-p channel dominant over the e-e one in the short- 
time regime (T < 100 fs) and for low energy carriers injected or optically pumped above the 
conduction band minimum but below the 2Eg threshold. This is due to the fact that one of 
the most important kinematic channels in the e-e scattering is the one that involves e-h pairs 
created across the energy gap, and those pairs are energetically forbidden when the carriers are 
pumped with energy that is 2Eg below the CBM. It is important to stress that on the basis 
of the relative balance of the e-e and e-p scattering channels presented in this work it is not 
possible be draw general conclusions for all kind of materials. If on one side I am convinced 
that, even if the CCA weakens the results obtained in the very-short time scale, the arguments 
based on the energy of the initial carrier to motivate a potential bigger importance of the e-p 
channel is reliable for any kind of semiconductor or insulator. Of course the same argument does 
no apply to metals. On the other hand silicon is not a system known for its strong e-p coupling. 
There are systems like bulk Diamond [TU], carbon nano-structures[10J or super-conductors [41) 
where the e-p scattering is enormously strong and this can strengthen the e-p channel over the 
e-e one. 
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